2017-18 Benthic Community Analysis

Purpose of This Document

Working with 2017-18 benthic invertebrate and fish diet data to calculate and plot:

Changes in scaper taxa community
NMS analysis
Diet NMS
Diet 1:1 and costello plot
Analysis of diet and benthic overlay
Effect of treatment on benthic and diet communities using MRPP or perMANOVA?

Benthic

Load Packages

# import data
bugs <- readxl::read_xlsx("./Data/2017-18 bugs.xlsx", sheet = "Cleaned")
ffg <- readxl::read_xlsx("./Data/2017-18 bugs.xlsx", sheet = "FFG")
sizes <- readxl::read_xlsx("./Data/2017-18 bugs.xlsx", sheet = "Size")
Diets <- readxl::read_xlsx("./Data/2018 Diets.xlsx")
Bentho <- read_xlsx("./Data/2018 Tiles Allison's Computer.xlsx", sheet = "Tiles")
Both_years<- read_excel("./Data/2018 Tiles Allison's Computer.xlsx", sheet = "2017 and 2018")

#Replace Ssp. in Taxon with Family
bugs$Taxon[bugs$Taxon == "Ssp."] <- bugs$Family[bugs$Taxon == "Ssp."]

#Setup ggplot theme
my_theme <- theme_bw(base_size = 18, base_family = "Microsoft Sans Serif") +
          theme(plot.title = element_text(hjust = 0.5),
                plot.caption = element_text(hjust = 0.5, face = "italic", color = "grey"),
                legend.title = element_text(colour = "black", size = 18, face = "bold"),
                legend.text = element_text(colour = "black", size = 12))
                
                
theme_set(my_theme)

Chlorophyll a

# Change variables to factors.  Re-order with CHUCK at the end
Bentho$Stream <- factor(Bentho$Stream, levels = c("LOON", "MCTE", "W-113", "W-100", "CHUCK", "W-122"))
Bentho$Treatment <- as.factor(Bentho$Treatment)
Bentho$Meter <- as.factor(Bentho$Meter)
Bentho$BenthoTotal <- as.numeric(Bentho$BenthoTotal)

# Drop W-122
Bentho <- Bentho %>% filter(Stream != "W-122")

# Calculate mean of the three tiles for each meter
bentho.mean <- aggregate(BenthoTotal ~ Stream + Treatment + Meter, mean, data = Bentho)

ggplot(data = bentho.mean) +
  geom_line(mapping = aes(x = Meter, y = BenthoTotal, color = Treatment, group = Treatment), size = 3) +
  ylab("Chlorophyll a (ug/cm^2)")+ scale_x_discrete(breaks = seq(0, 90, by = 15)) +
  scale_color_viridis_d(option = "E") +
  facet_wrap("Stream")

# Change variables to factors.  Re-order with CHUCK at the end
Both_years$Stream <- factor(Both_years$Stream, levels = c("LOON", "MCTE", "W-113", "W-100", "CHUCK", "W-122"))
Both_years$Year <- factor(Both_years$Year, levels = c("2018", "2017"))
Both_years$Year.Treatment <- as.factor(Both_years$Year.Treatment)
Both_years$Treatment <- as.factor(Both_years$Treatment)
Both_years$Meter <- as.factor(Both_years$Meter)
Both_years$BenthoTotal <- as.numeric(Both_years$BenthoTotal)

# Drop W-122.
Both_years <- Both_years %>% filter(Stream != "W-122")

# Calculate mean of the three tiles for each meter
bentho_both.mean <- aggregate(BenthoTotal ~ Stream + Treatment + Year.Treatment + Meter + Year, mean, data = Both_years)
bentho.reach <- aggregate(BenthoTotal ~ Stream + Treatment + Year.Treatment + Year, mean, data = Both_years)
write_csv(bentho.reach, "./Data/bentho.reach.csv")

ggplot(data = bentho_both.mean) +
  geom_line(mapping = aes(x = Meter, y = BenthoTotal, color = Treatment, 
  linetype = Year, group = Year.Treatment), size = 3) +
  ylab("Chlorophyll a (ug/cm^2)") + 
  scale_x_discrete(breaks = seq(0, 90, by = 15)) +
  scale_color_viridis_d(option = "E") +
  facet_wrap("Stream") +
  theme(legend.position = c(.84, .27), legend.key.width = unit(5, "cm"))

Setup

In order to calculate density we must take into account that our count values are from a subsample of three pooled samples.

#Calculations...
bugs$Density <- bugs$Count / bugs$PercentSub / .09  #Calculate total in pooled sample

bugs$CollDate <- as.Date(bugs$CollDate, format = "%m/%d/%y") %>% year() %>% as.factor()

bugs.agg <- bugs %>%
  group_by(CollDate, Stream, Treatment, Taxon) %>%
  summarise_at(vars(Density), funs(sum)) %>% ungroup

bugs.agg$Density <- bugs.agg$Density / 3    #Divide by number of samples pooled

We divide Count by percent subsampled (actually a fraction) to get a count for the total sample taken. We then divide by .09 which is the area of the surber sampler (in m\(^2\)). We then aggregate and divide by three (There were three samples taken per reach).

# Spread and then gather to fill missing Taxons with 0's
bugs.agg <- spread(bugs.agg, key = "Taxon", value = "Density") %>% #bugs.agg is benthic samples
  mutate_if(is.numeric , replace_na, replace = 0) %>% 
  gather(key = "Taxon", value = "Density", 4:ncol(.))

# Add FFG and size classes
bugs.agg <- merge(bugs.agg, sizes) %>% merge(., ffg)

# Create list of CollDate, Stream and Taxon
bugs.reach.diff <- bugs.agg %>% select(-c("Density", "CollDate")) %>% unique()

# Calculate Differences
bugs.reach.diff$Diff <- bugs.agg$Density[bugs.agg$CollDate == "2018"] -     
  bugs.agg$Density[bugs.agg$CollDate == "2017"]

# Filter for only 2018 bugs
bugs18 <- filter(bugs.agg, CollDate == "2018")

A quick taste of the finished data table product:

# Datatable with filters and data ranges
datatable(bugs.agg, rownames = FALSE, class = "hover", filter = "top", options = list(pageLength = 10, scrollX = T)) %>% formatRound("Density", digits = 2, interval = 3, 
                                          mark = ",", dec.mark = getOption("OutDec"))

Taxa

# Plot density differences
ggplot(data = bugs.reach.diff, aes(x = Taxon, y = Diff, fill = Treatment)) + 
  geom_col(position = "dodge") +
  facet_wrap(c("Stream")) +
  ggtitle("Differences in Taxa Abundances Between Reaches Pre and Post Gap Years") +
  scale_fill_viridis_d(option = "E") +
  coord_flip()

Plot of density differences between 2017 and 2018 for both reaches. As seen, there are no ubiquitous taxa that consistently increase or decrease across streams.

FFG’s

#Aggregate to get density differences of FFG's
bugs.ffg <- bugs.reach.diff %>% 
  group_by(FFG, Stream, Treatment) %>%
  summarise_at(vars(Diff), funs(sum)) %>% ungroup

#Add total SC, this throws off counts (scrapers doubly reprented)
ffg.SC <- filter(bugs.ffg, FFG %in% c("SCi", "SCe"))
ffg.SC$FFG <- "SC"
ffg.SC <- ffg.SC %>% group_by(FFG, Stream, Treatment) %>%
  summarise_at(vars(Diff), funs(sum)) %>% ungroup
ffg.SC <- rbind(bugs.ffg, ffg.SC)

ggplot(data = ffg.SC, aes(x = FFG, y = Diff, color = Treatment, shape = Stream)) + 
  geom_jitter(width = .2, height = .2, size = 4) +
  labs(title = "FFG Differences Between Years (Treatment and Control Reaches)", 
       caption = "CF = Collector Filterer, CG = Collector Gatherer,  SC = Scrapers (a composite of SCe and SCi), \n         SCe = Scrapers that are edible, SCi = Scrapers that are inedible, P = Predators") +
  scale_color_viridis_d(option = "E") +
  coord_flip()

Here we plot density differences of FFG’s.

Notice the only group with consistently elevated differences in the post gap year (2018) is Scrapers.

bugs18 %>%
  group_by(Stream, Treatment, FFG) %>%
  summarise_at(vars(Density), funs(sum)) %>% 
  ungroup() %>%
  
  ggplot(aes(x = FFG, y = Density, color = Treatment, shape = Stream)) + 
    geom_jitter(width = .2, height = .2, size = 4) + 
    ggtitle("2018 FFG Densities (Control & Treatment)")  +
    scale_color_viridis_d(option = "E") +
    coord_flip()

Here we can see that Collector Gatherers are the most abundant taxa, but the aggregate of the two scraper categories (SCe and SCi) is also up there.

When comparing between the reaches we see elevated abundances of both scraper groups (except for SCi in CHUCK). Collector Gatherers are also consistently elevated (except in CHUCK), which seems contradictory to the between year comparison (Control reach has always had elevated abundances of CG’s?)

Will these functional feeding groups be important in diets?
Does the relative change in abundance of a taxa group change fish selection (i.e. now that there are more scrapers in the treatment reach are the fish going to town on Scrapers?)
TBD. (see other.md file in the Diets project)

bugs.reach.diff %>% filter(FFG == "SCe" | FFG == "SCi") %>% 
  ggplot(aes(x = Taxon, y = Diff, color = Treatment, shape = Stream)) + 
    geom_jitter(width = .2, height = .2, size = 4) +
    ggtitle("Diff. Between Reaches of Scraper Taxa (2017 & 2018") +
    scale_color_viridis_d(option = "E") +
    coord_flip()

Here we see that the Scraper taxa with the greatest change are:

Micrasema
Juga
Glossosoma
Drunella
Heptageniidae taxa

bugs18 %>% filter(FFG == "SCe" | FFG == "SCi") %>%
  ggplot(aes(x = Taxon, y = Density, color = Treatment, shape = Stream)) +
    geom_jitter(width = .2, height = .2, size = 4) +
    ggtitle("Density of Scraper Taxa in 2018") +
    scale_color_viridis_d(option = "E") +
    coord_flip()

bugs.reach.diff %>% 
  group_by(Size, Stream, Treatment) %>%
  summarise_at(vars(Diff), funs(sum)) %>%
  ungroup() %>%
    ggplot(aes(x = Size, y = Diff, color = Treatment, shape = Stream)) + 
      geom_jitter( width = .2, height = .2, size = 4) +
      ggtitle("Differences in Size Classes") +
      scale_color_viridis_d(option = "E")

Two size classes are plotted, small (S) and large (L) based on nothing, just observation… Hopefully I will soon incorporate some information from the Poff database about final larval instar size and have three size classes: small, medium, and large.

NMS

Who knows what I’m actually doing, all of this ordination stuff is currently pre BOT 570. I will definitely update all this once I know how to handle singleton taxa, loads of zero’s, etc.

# Spread data to get a taxon matrix, create factor of CollDate and Treatment
bugs.mds <- bugs.agg %>% 
  select(-c("FFG", "Size")) %>%
  group_by(Stream, Treatment, CollDate, Taxon) %>%  
  spread(key = "Taxon", value = "Density") %>%
  mutate_if(is.numeric , replace_na, replace = 0) %>% ungroup() %>%
  mutate(YearTreat = interaction(CollDate, Treatment))

# Get grouping variables and taxon matrix
Enviro <- bugs.mds %>%
  select(Stream, Treatment, CollDate, YearTreat) 

Density <- bugs.mds %>% 
  select("Ameletus":"Zapada")

grouping <- select(Enviro, Stream, Treatment, CollDate, YearTreat)

#Get MDS for invert density
sol <- metaMDS(Density, distance = "bray", k = 2, trymax = 100)

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1726121 
## Run 1 stress 0.1726121 
## ... New best solution
## ... Procrustes: rmse 1.088343e-05  max resid 2.412069e-05 
## ... Similar to previous best
## Run 2 stress 0.1726121 
## ... Procrustes: rmse 5.145144e-06  max resid 1.288155e-05 
## ... Similar to previous best
## Run 3 stress 0.1726121 
## ... Procrustes: rmse 2.978185e-05  max resid 7.206061e-05 
## ... Similar to previous best
## Run 4 stress 0.2287457 
## Run 5 stress 0.1726122 
## ... Procrustes: rmse 0.0001131157  max resid 0.0002641682 
## ... Similar to previous best
## Run 6 stress 0.1726121 
## ... Procrustes: rmse 5.036715e-05  max resid 0.0001125872 
## ... Similar to previous best
## Run 7 stress 0.1726121 
## ... New best solution
## ... Procrustes: rmse 2.371582e-06  max resid 4.925835e-06 
## ... Similar to previous best
## Run 8 stress 0.1726122 
## ... Procrustes: rmse 0.0001348066  max resid 0.0003069075 
## ... Similar to previous best
## Run 9 stress 0.1726122 
## ... Procrustes: rmse 0.0001038851  max resid 0.0002312197 
## ... Similar to previous best
## Run 10 stress 0.1726121 
## ... Procrustes: rmse 5.222521e-06  max resid 1.278906e-05 
## ... Similar to previous best
## Run 11 stress 0.1726121 
## ... Procrustes: rmse 3.285539e-05  max resid 8.987912e-05 
## ... Similar to previous best
## Run 12 stress 0.1726121 
## ... Procrustes: rmse 4.742137e-05  max resid 0.0001067743 
## ... Similar to previous best
## Run 13 stress 0.1726121 
## ... Procrustes: rmse 1.159735e-05  max resid 2.82002e-05 
## ... Similar to previous best
## Run 14 stress 0.2437675 
## Run 15 stress 0.1726121 
## ... Procrustes: rmse 6.751397e-05  max resid 0.0001508694 
## ... Similar to previous best
## Run 16 stress 0.1726121 
## ... Procrustes: rmse 1.473293e-05  max resid 3.183357e-05 
## ... Similar to previous best
## Run 17 stress 0.1726121 
## ... Procrustes: rmse 2.726293e-05  max resid 6.344963e-05 
## ... Similar to previous best
## Run 18 stress 0.2488585 
## Run 19 stress 0.1726121 
## ... Procrustes: rmse 1.263578e-05  max resid 2.885207e-05 
## ... Similar to previous best
## Run 20 stress 0.1738211 
## *** Solution reached

This is a pretty good stress value, who knows how that will change when I start doing this right… Can see the transfrmations applied at the top of the output, this was done automatically, I wouldn’t know what to do. Used the Bray-Curtis distance metric, max tries is set to 100 with dimensions equal to 2.

#set up NMDS with dimensions of sol and env factors from "Enviro". 
NMDS <- data.frame(x = sol$points[, 1], y = sol$points[ ,2], 
                   Stream = select(grouping, Stream), 
                   Treatment = select(grouping, Treatment), 
                   CollDate = select(grouping, CollDate),
                   YearTreat = select(grouping, YearTreat))

#Get mean x,y values 
NMDS.mean <- select(NMDS, x, y) %>% aggregate(grouping, mean)

#Load ellipse for NMDS
plot.new()
ord <- ordiellipse(sol, NMDS$YearTreat, display = "sites", kind = "sd", conf = 0.95, label = TRUE)

dev.off()

## null device 
##           1

The “display” argument can be set to “site” or “species” depending on what you want to group. “kind” can be standard deviation or standard error, not sure which to use here.

# Fit the ellipse function to actual data
df_ell <- data.frame()
for(g in NMDS$YearTreat){
  df_ell <- rbind(df_ell, cbind(as.data.frame(with(NMDS[NMDS$YearTreat == g,],
                  vegan:::veganCovEllipse(ord[[g]]$cov, ord[[g]]$center, ord[[g]]$scale)))
                                ,YearTreat = g))
}

not going to pretend like I know how this works, got it from https://stackoverflow.com/questions/13794419/plotting-ordiellipse-function-from-vegan-package-onto-nmds-plot-created-in-ggplo but I do know that changing the column selected from NMDS changes which variable is used for grouping.

ggplot(data = NMDS, aes(x, y, color = YearTreat)) +
  geom_path(data = df_ell, aes(x = NMDS1, y = NMDS2)) +
  annotate("text", x = (NMDS$x), y = (NMDS$y) + .04, label = NMDS$Stream, size = 2) +
  geom_point(aes(shape = Stream)) +
  ggtitle("NMDS of Benthic Invertebrate Community") +
  scale_color_viridis_d()

Actually maybe cool result, all of the treatment reaches now plot closer post-gap year. The controls show a similar pattern though maybe a little less pronounced. Everything also shifted left between the two years, although this isn’t because of the treatment, just annual changes. Maybe the gap amplified whatever effect the year had, extra, extra light or temperature maybe.

Summary Info

bugs.agg %>% mutate(CollDate.Treatment = interaction(Treatment, CollDate)) %>%
group_by(Stream, CollDate.Treatment) %>%
  summarise_at(vars(Density), funs(sum)) %>%
  ggplot(aes(x = Stream, y = Density, fill = CollDate.Treatment)) + 
  geom_col(position = "dodge") +
  ggtitle("Total Inverebrate Density (m^2)") +
  scale_fill_viridis_d(option = "E")

  #datatable(rownames = FALSE, class = "hover", filter = "top", 
            #options = list(pageLength = 10, scrollX=T), ) %>% 
            #formatRound("Density", digits = 2, interval = 3, mark = ",", dec.mark = getOption("OutDec"))

group_by(bugs.reach.diff, Stream, Treatment) %>%
  summarise_at(vars(Diff), funs(sum)) %>% 
  arrange(Treatment) %>% 
  ggplot(aes(x = Stream, y = Diff, fill = Treatment)) + 
  geom_col(position = "dodge") +
  ggtitle("Differences in Density (T-C)") +
  scale_fill_viridis_d(option = "E")

  #datatable(rownames = FALSE, class = "hover", filter = "top", 
            #options = list(pageLength = 10, scrollX=T), ) %>% 
            #formatRound("Diff", digits = 2, interval = 3, mark = ",", dec.mark = getOption("OutDec"))

bugs.agg %>% mutate(CollDate.Treatment = interaction(Treatment, CollDate)) %>%
  group_by(Stream, CollDate.Treatment, FFG) %>%
  summarise_at(vars(Density), funs(sum)) %>% write_csv("./Data/FFG.Densities.csv")
  
  
  #ggplot(aes(x = FFG, y = Density, fill = CollDate.Treatment)) + 
  #geom_col(position = "dodge") +
  #facet_wrap("Stream") +
  #ggtitle("Differences in Density By FFG (T-C)") +
  #theme_bw() +
     #theme(plot.title = element_text(size = 40, face = "bold")) +
  #scale_fill_viridis_d(option = "E") 
  
  #arrange(CollDate) %>% write_csv("./Data/ffg.diffs.csv")
  
  #datatable(rownames = FALSE, class = "hover", filter = "top", options = list(pageLength = 10, #scrollX=T), ) %>% formatRound("Diff", digits = 2, interval = 3, 
                                         #mark = ",", dec.mark = getOption("OutDec"))

Diets

# Fill missing family's with order
Diets$Family[is.na(Diets$Family)] <- Diets$Order[is.na(Diets$Family)]

# Merge diets, sizes and FFG's
Diets <- merge(Diets, sizes, by.x = "Family", by.y = "Taxon") %>% 
  merge(., ffg, by.x = "Family", by.y = "Taxon")

# Aggregate by rep
Diets.by.fish <- Diets %>% 
  group_by(Family, Stream, Treatment, Origin, FFG, Size, Rep) %>%
  summarise_at("Count", funs(sum)) %>% ungroup

# Aggregate by reach
Diets.reach <- Diets %>% 
  group_by(Family, Stream, Treatment, Origin, FFG, Size) %>%
  summarise_at("Count", funs(sum)) %>% ungroup

bugs.family <- bugs %>% 
  filter(CollDate == "2018") %>%
  select(-Taxon) %>% 
  group_by(Stream, Treatment, Family, CollDate) %>%
  summarise_at(vars(Density), funs(sum)) %>% ungroup

bugndiet <- merge(Diets.reach, bugs.family, all.x = TRUE, all.y = TRUE) %>% 
  filter(Origin == "A") %>%
  mutate_if(is.numeric, replace_na, replace = 0)

Diets.reach$Stream.Treat <- interaction(Diets.reach$Treatment, Diets.reach$Stream)

Diets.reach %>%
  group_by(Stream.Treat, Origin) %>%
  summarise_at(vars(Count), funs(sum)) %>% ungroup() %>%

ggplot(aes(x = Stream.Treat, y = Count, fill = Origin)) +
  geom_col() +
  labs(title = "Total Aquatics and Terrestricals by Reach") +
    theme(legend.position = c(.84, .8), legend.key.width = unit(5, "cm")) +
    scale_fill_viridis_d(option = "E", end = .8)

?scale_fill_viridis_d

By Family

# Plot proportion of a taxon in the benthic community versus proporion of a taxon in the diet community

bugndiet %>% 
  group_by(Stream, Treatment, Family) %>% 
  summarise_if(is.numeric, funs(sum)) %>% 
  mutate(frac.diet = Count / sum(Count)) %>%
  mutate(frac.benthic = Density / sum(Density)) %>% ungroup() %>%
  
  ggplot(aes(x = frac.benthic, y = frac.diet, color = Treatment)) +
    geom_point() +
    geom_abline(slope = 1, intercept = 0) +
    coord_cartesian(ylim = c(0, .8), xlim = c(0, .8)) +
    geom_text(aes(x = frac.benthic, y = frac.diet, label = Family), size = 3, nudge_y = .02) +
    ggtitle("1-1 Prop. of Taxa in Benthic Comm. Vs. Prop. in Diets") +
    facet_wrap("Stream") +
    scale_color_viridis_d(option = "E") +
    theme(legend.position = c(.84, .27), legend.key.width = unit(5, "cm"))

Plot of proportional abundance of an individual taxon in the benthic community versus the proportional abundance of the aggregate diets for that reach.

Here we see that Chironomidae typically make up the greatest abundance of both the fish diets and the benthic community, with LOON being the exception (Brachycentridae are equally represented in the diets of control reach and over represented in the treatment). In MCTE, Juga compose a large part of the diets relative to their abundance due to one outlier fish.

There are no overarching differences in how a taxa maps out given the treatment, but the difference in Brachycentridae in LOON between treatment and control is interesting.

By FFG

# Plot proportion of a FFG in the benthic community versus proporion of a FFG in the diet community
bugndiet %>% 
  group_by(Stream, Treatment, FFG) %>% 
  summarise_if(is.numeric, funs(sum)) %>% 
  mutate(frac.diet = Count / sum(Count)) %>%
  mutate(frac.benthic = Density / sum(Density)) %>%
  mutate_if(is.character, replace_na, replace = "T") %>%

  ggplot(aes(x = frac.benthic, y = frac.diet, color = Treatment)) +
    geom_point() +
    geom_abline(slope = 1, intercept = 0) +
    coord_cartesian(ylim = c(0, 1), xlim = c(0, 1)) +
    geom_text(aes(x = frac.benthic, y = frac.diet, label = FFG), size = 4, nudge_y = .02) +
    facet_wrap("Stream") +
    ggtitle("1-1 Prop. of FFG's in Benthic Comm. Vs. Prop. in Diets") +
    scale_color_viridis_d(option = "E") +
    theme(legend.position = c(.84, .27), legend.key.width = unit(5, "cm"))

Plot of proportional abundance of an individual FFG in the benthic community versus the proportional abundance of the aggregate diets for that reach.

We see that Collector Gatherers and Shredders are the most abundant both in diets and in the benthic community. There aren’t any over arching trends in how FFG’s fall out by treatment.

By Size Class

# Plot proportion of a size class in the benthic community versus proporion of a FFG in the diet community
bugndiet %>% 
  group_by(Stream, Treatment, Size) %>% 
  summarise_if(is.numeric, funs(sum)) %>% 
  mutate(frac.diet = Count / sum(Count)) %>%
  mutate(frac.benthic = Density / sum(Density)) %>%
  mutate_if(is.character, replace_na, replace = "T") %>%

  ggplot(aes(x = frac.benthic, y = frac.diet, color = Treatment)) +
    geom_point() +
    geom_text(aes(x = frac.benthic, y = frac.diet, label = Size), size = 4, nudge_y = .02) +
    geom_abline(slope = 1, intercept = 0) +
    coord_cartesian(ylim = c(0, 1), xlim = c(0, 1)) +
    facet_wrap("Stream") +
    ggtitle("1-1 Prop. of Size Classes in Benthic Comm. Vs. Prop. in Diets") +
    scale_color_viridis_d(option = "E") +
    theme(legend.position = c(.84, .27), legend.key.width = unit(5, "cm"))

Plot of proportional abundance of an individual size class in the benthic community versus the proportional abundance of the aggregate diets for that reach.

Costello method plot

Diets.fish <- Diets.by.fish %>% 
  group_by(Stream, Treatment, Rep) %>%
  mutate(frac.diet = Count / sum(Count))

Diets.fish <- Diets.fish %>% 
  group_by(Stream, Treatment, Family) %>%
  mutate(frac.fish = length(Family[Count > 0])) %>%
  group_by(Stream, Treatment) %>%
  mutate(frac.fish = frac.fish / max(Rep)) %>%
  group_by(Stream, Treatment, Family)

Diets.fish <- Diets.fish %>%
  summarise_at(vars(frac.fish, frac.diet), funs(mean))

ggplot(data = Diets.fish) +
  geom_point(mapping = aes(x = frac.fish, y = frac.diet, color = Treatment)) +
  geom_text(aes(x = frac.fish, y = frac.diet, label = Family), size = 3, nudge_y = .02) +
  geom_abline(slope = 1, intercept = 0) +
  coord_cartesian(ylim = c(0, 1), xlim = c(0, 1)) +
  facet_wrap("Stream") +
  ggtitle("Costello Plot of Prop. Fish Occurance Vs. Prop. Diet Community") +
  scale_color_viridis_d(option = "E") +
  theme(legend.position = c(.84, .27), legend.key.width = unit(5, "cm"))

The Costello method plots % occurance of a taxon in fish versus % of aggregate diet. With the adjusted method the aggregate diet is only of fish that had the taxon present (ignoring zero values). We see that in MCTE and in W-100 there seem to be just a couple fish that specifically target Juga and ELmidae respectively. In MCTE, the taxa from treatment diets seem to consistently constitute a smaller portion of the diet than would be expected given how many fish they occur in.

Diet MDS